# 三角分解法
import numpy as np


def lufact(A, B):
    n = B.size
    Y = np.zeros([n, 1], dtype=np.double)
    R = np.arange(0, n)
    epslion = np.finfo(np.float32).eps

    for p in range(n):
        max_index = np.argmax(np.abs(A[p:n, p]))
        max_row = max_index + p
        A[[p, max_row], :] = A[[max_row, p], :]
        # R[p], R[max_row] = R[max_row], R[p]
        R[[p, max_row]] = R[[max_row, p]]

        if abs(A[p, p]) < epslion:
            print('A is singular. No unique solution')
            break

        for k in range(p + 1, n):
            mult = A[k, p] / A[p, p]
            A[k, p] = mult
            A[k, p + 1:n] = A[k, p + 1:n] - mult * A[p, p + 1:n]

    Y[0] = B[R[0]]
    for k in range(1, n):
        Y[k] = B[R[k]] - A[k, 0:k] @ Y[0:k]

    return A, Y


def backsub(A, B):
    n = B.size

    X = np.zeros([n, 1], dtype=np.double)
    X[n - 1] = B[n - 1] / A[n - 1][n - 1]
    for i in range(n - 2, -1, -1):
        X[i] = (B[i] - A[i, i + 1:] @ X[i + 1:]) / A[i][i]

    return X


def out(x):
    print(x)


def main():
    A = np.array([[1 / 4, 1 / 5, 1 / 6],
                  [1 / 3, 1 / 4, 1 / 5],
                  [1 / 2, 1, 2]], dtype=np.float64)
    # 解
    B = np.array([9, 8, 8], dtype=np.float64).reshape(3, 1)

    A, B = lufact(A, B)
    out(A)

    X = backsub(A, B)
    out(X)


if __name__ == '__main__':
    main()
